library(funFEM)
library(reshape2)
library(ggplot2)
library(gridExtra)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(mclust)
library(stringr)
library(cluster)
library(clusterSim)
library(ggmap)
Cet ensemble de données contient des données provenant du système de partage de vélos de Paris, appelé Vélib. Les données sont des profils de chargement des stations de vélos sur une semaine. Les données ont été collectées toutes les heures pendant la période du dimanche 1er septembre au dimanche 7 septembre 2014.
On charge ici les données Velib disponibles dans la
librairie funFEM.
library(funFEM)
data(velib)
Ces données contiennent
data : les profils de chargement (nb de vélos disponibles / nb de points d’attache) des 1189 stations à 181 points de temps.
position : la longitude et la latitude des 1189 stations de vélos.
dates : les dates de téléchargement.
bonus : indique si la station est sur une colline (bonus = 1).
names : les noms des stations.
Afin d’avoir des semaines complètes, nous allons supprimer les 13 premières colonnes.
Velibdata <- velib$data[, -c(1:13)]
colnames(Velibdata) <- velib$dates[-c(1:13)]
On controle le nombre de jours et d’heures dans le jeu de données
day <- str_sub(colnames(Velibdata), 1, 3)
day <- factor(day, levels = c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim"))
table(day)
day
Lun Mar Mer Jeu Ven Sam Dim
24 24 24 24 24 24 24
hour <- str_sub(colnames(Velibdata), 5, 6)
table(hour)
hour
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
timeTick <- 24 * (0:6)
On stocke les informations sur les stations dans le dataFrame suivant :
station <- data.frame(nom = velib$names, velib$position, colline = velib$bonus)
On ajoute des informations complémentaires sur les stations :
# Les gares
gare <- rep(0, nrow(station))
gare[which(str_detect(station$nom, "GARE") == T)] <- 1
# Lieux touristiques (autour de Tour eiffel, Panthéon, Arc de Triomphe, Champs
# Elysées, Louvre, Notre Dame, Invalides, Trocadéro, Sacré Coeur, Montmartre)
tourisme <- rep(0, nrow(station))
TI <- c(209, 504, 897, 102, 518, 357, 256, 301, 992, 627, 579, 241, 646, 13, 283,
941, 976, 1189, 22, 633, 853, 973, 762, 195, 71, 131, 912, 302, 303, 1068, 774,
384, 307, 326, 436, 789, 250, 755, 461, 439, 918, 1065, 99, 151, 313, 871)
tourisme[TI] <- 1
# Les mairies
mairie <- rep(0, nrow(station))
mairie[which(str_detect(station$nom, "MAIRIE") == T)] <- 1
mairie[which(str_detect(station$nom, "HOTEL") == T)] <- 1
# Bois / Parc / Square
parc <- rep(0, nrow(station))
parc[which(str_detect(station$nom, "BOIS DE") == T)] <- 1
parc[which(str_detect(station$nom, "PARC") == T)] <- 1
parc[which(str_detect(station$nom, "SQUARE") == T)] <- 1
parc[which(str_detect(station$nom, "CANAL") == T)] <- 1
# Porte
porte <- rep(0, nrow(station))
porte[which(str_detect(station$nom, "PORTE") == T)] <- 1
station <- data.frame(station, gare = gare, tourisme = tourisme, mairie = mairie,
parc = parc, porte = porte)
Fonction pour tracer les profils (charge de chaque station / loading) en fonction des jours et heures pour les stations choisies.
plotprofils <- function(Velibdata, station, numstation, plotsem = TRUE) {
# Velibdata = les données de velib initiales station = ens. des
# informations sur les stations numstation = vecteur contenant le numéro de
# ligne des stations dont on souhaite tracer les profils
library(reshape2)
library(ggplot2)
Dataaux <- data.frame(id.s = station$nom, Velibdata)
Aux <- melt(Dataaux[numstation, ], id = c("id.s"))
Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
5, 6)))
Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
"Ven"))
if (plotsem == TRUE) {
g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
ncol = 2) + xlab("Time") + ylab("Loading") + scale_x_continuous(breaks = seq(0,
24, 4)) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
labs(col = "Classif.") + scale_x_continuous(breaks = c(0, 7, 9, 12, 14,
20, 24))
} else {
v <- rep(0, nrow(Aux))
for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
levels(Aux$day)[j])] + (24 * (j - 1))
Aux <- data.frame(Aux, v = v)
rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
ymax = rep(1, 7), color = rainbow(n = 7))
g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
ylab("Loading") + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
}
return(g)
}
Exemple d’utilisation :
plotprofils(Velibdata, station, numstation = c(813, 655, 468))
plotprofils(Velibdata, station, numstation = c(813, 655, 468), plotsem = FALSE)
Se donnant une classification, fonction pour tracer les profils moyens par classe :
plotprofilsmoy <- function(Velibdata, clustering, plotsem = TRUE) {
library(reshape2)
library(ggplot2)
Dataaux <- matrix(0, nrow = max(clustering), ncol = ncol(Velibdata))
for (k in 1:max(clustering)) {
Dataaux[k, ] <- apply(Velibdata[which(clustering == k), ], 2, mean)
}
colnames(Dataaux) <- colnames(Velibdata)
Dataaux <- data.frame(id.s = as.factor(1:max(clustering)), Dataaux)
Aux <- melt(Dataaux, id = c("id.s"))
Aux <- data.frame(Aux, day = str_sub(Aux$variable, 1, 3), hour = as.numeric(str_sub(Aux$variable,
5, 6)))
Aux$day <- factor(Aux$day, levels = c("Lun", "Mer", "Mar", "Sam", "Jeu", "Dim",
"Ven"))
if (plotsem == TRUE) {
g <- ggplot(Aux, aes(x = hour, y = value, col = id.s)) + geom_line() + facet_wrap(~day,
ncol = 2) + xlab("Time") + ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") +
theme(axis.title.x = element_blank()) + labs(col = "Classif.") + scale_x_continuous(breaks = c(0,
7, 9, 12, 14, 20, 24))
} else {
v <- rep(0, nrow(Aux))
for (j in 1:length(levels(Aux$day))) v[which(Aux$day == levels(Aux$day)[j])] <- Aux$hour[which(Aux$day ==
levels(Aux$day)[j])] + (24 * (j - 1))
Aux <- data.frame(Aux, v = v)
rect <- data.frame(xmin = timeTick, xmax = timeTick + 23, ymin = rep(0, 7),
ymax = rep(1, 7), color = rainbow(n = 7))
g <- ggplot(data = Aux) + geom_line(data = Aux, aes(x = v, y = value, col = id.s)) +
geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
fill = color), show.legend = FALSE, alpha = 0.1) + xlab("Time") +
ylab("Loading") + ylim(0, 1) + theme(legend.position = "bottom") + theme(axis.title.x = element_blank()) +
labs(col = "Classif.") + scale_x_continuous(breaks = rep(seq(0, 21, 3),
7) + rep(24 * (0:6), each = 4), labels = rep(seq(0, 21, 3), 7))
}
return(g)
}
Exemple d’utilisation avec les stations sur colline par rapport aux autres :
plotprofilsmoy(Velibdata, clustering = station$colline + 1, plotsem = FALSE)
Fonction auxiliaire pour comparer une classification avec les informations auxiliaires disponibles pour les stations
stationcaract <- function(station, clustering) {
library(dplyr)
# Dataaux<-data.frame(id.s=c(1:nrow(Data)),station) aux <- cbind(Dataaux,
# clust)
aux <- cbind(station, clustering)
aux.long <- melt(data.frame(lapply(aux[, -c(2:3)], as.character)), stringsAsFactors = FALSE,
id = c("nom", "clustering"), factorsAsStrings = T)
# Effectifs
aux.long.q <- aux.long %>%
group_by(clustering, variable, value) %>%
mutate(count = n_distinct(nom)) %>%
distinct(clustering, variable, value, count)
# avec fréquences
aux.long.p <- aux.long.q %>%
group_by(clustering, variable) %>%
mutate(perc = count/sum(count)) %>%
arrange(clustering)
gaux <- ggplot(data = aux.long.p, aes(x = clustering, y = perc, fill = value)) +
geom_bar(stat = "identity") + facet_wrap(~variable)
return(gaux)
}
On reprend ici rapidement quelques éléments d’analyse descriptive vus précédemment dans l’UF.
library(corrplot)
timeTickAux <- c(timeTick, 168)
for (k in 1:7) corrplot(cor(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k + 1])]),
method = "ellipse")
velibmeanday <- matrix(0, nrow = nrow(Velibdata), ncol = 7)
for (k in 1:7) {
velibmeanday[, k] <- apply(Velibdata[, (timeTickAux[k] + 1):(timeTickAux[k +
1])], 1, mean)
}
colnames(velibmeanday) <- c("Lun", "Mar", "Mer", "Jeu", "Ven", "Sam", "Dim")
ggplot(melt(velibmeanday), aes(x = Var2, y = value)) + geom_boxplot() + xlab("") +
ylab("Loading")
velibHour <- data.frame(value = colMeans(Velibdata), day = day, hour = as.numeric(hour))
ggplot(velibHour, aes(x = hour, y = value, col = day)) + geom_line() + geom_point() +
ylab("Loading") + ggtitle("Hourly loading, averaged over all stations")
library(ggmap)
# Average loading for each station
load <- rowMeans(Velibdata)
qmplot(longitude, latitude, data = station, color = load) + scale_colour_gradient(high = "red",
low = "yellow")
# scale_colour_gradient(high='blue',low='white')
# Loading moyen à l'heure t (le premier temps est 0h)
t <- 0
loadheure <- rowMeans(Velibdata[, seq(from = (t + 1), by = 24, length = 7)])
qmplot(longitude, latitude, data = station, color = loadheure) + scale_colour_gradient(high = "red",
low = "yellow")
qmplot(longitude, latitude, data = station, color = as.factor(station$colline))
qmplot(longitude, latitude, data = station, color = as.factor(station$tourisme))
Question : Faites une ACP des stations et analysez les résultats.
respca <- PCA(.......)
fviz_eig(....)
fviz_pca_var(......)
fviz_pca_ind(........)
Dans cette partie, nous allons utiliser les coordonnées de l’ACP comme matrice de données pour classer les stations.
velibacp <- respca$ind$coord[, 1:7]
Question : A l’aide d’une procédure des Kmeans, déterminez une classification des stations.
# A COMPLETER
Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue.
# A COMPLETER
Question : A l’aide d’une procédure hiérarchique, déterminez une classification des stations.
# A COMPLETER
Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par CAH. Comparez cette classification avec celle obtenue par les Kmeans.
# A COMPLETER
Question : A l’aide des mélanges gaussiens, déterminez une classification des stations.
# A COMPLETER
Question : A l’aide des fonctions auxiliaires données dans la première partie, étudiez la classification obtenue par modèles de mélange. Comparez cette classification avec celles obtenues précédemment.
# A COMPLETER
Dans cette partie, on va s’intéresser à un clustering sur un jeu de données “résumé” au sens suivant : on fait la moyenne des loadings selon les tranches horaires suivantes : 0-7h, 8h-20h et 21h-23h pour chaque jour et chaque station.
VelibMean <- matrix(0, nrow = nrow(Velibdata), ncol = 3 * 7)
for (jour in 1:7) {
VelibMean[, 3 * (jour - 1) + 1] <- apply(Velibdata[, 24 * (jour - 1) + c(1:8)],
1, mean)
VelibMean[, 3 * (jour - 1) + 2] <- apply(Velibdata[, 24 * (jour - 1) + c(9:21)],
1, mean)
VelibMean[, 3 * (jour - 1) + 3] <- apply(Velibdata[, 24 * (jour - 1) + c(22:24)],
1, mean)
}
colnames(VelibMean) <- paste(rep(unique(day), each = 3), "-", rep(c("0-7h", "8h-20h",
"21-23h"), 7), sep = "")
Question : Etudiez les corrélations entre variables
du jeu de données VelibMean.
# A COMPLETER
Question : Faites une ACP et commentez.
# A COMPLETER
Question : Faites une classification des stations à
partir du jeu de données VelibMean à l’aide des Kmeans.
Etudiez la classification obtenue.
# A COMPLETER
Question : Faites une classification des stations à
partir du jeu de données VelibMean à l’aide des mélanges
gaussiens. Etudiez la classification obtenue.
# A COMPLETER
Les données de Vélib sont des données fonctionnelles (on étudie un
phénomène qui évolue au cours du temps). Il existe des méthodes de
clustering dédiées aux données fonctionnelles, comme par exemple le
package funFEM.
Pour utiliser cette méthode, on commence par considérer la décompostion dans la base de Fourier des courbes.
basis <- create.fourier.basis(c(0, 167), nbasis = 25)
fdobj <- smooth.basis(1:168, t(Velibdata), basis)$fd
La méthode funFEM est une procédure de clustering basée
sur des mélanges fonctionnels. Comme pour les mélanges gaussiens, il y a
différentes formes disponibles selon les contraintes imposés dans le
modèle (le détail est hors programme !). Cette méthode est détaillée
dans l’article de Bouveyron et al. (2015).
Question : Executez les commandes suivantes et exploitez les résultats. On considère une classification en 10 classes comme dans Bouveyron et al. (2015).
resfunFEM <- funFEM(fdobj, K = 10, model = "AkjBk", init = "kmeans", lambda = 0,
disp = TRUE)
df <- data.frame(proba = apply(resfunFEM$P, 1, max), classif = as.factor(apply(resfunFEM$P,
1, which.max)))
ggplot(df, aes(x = classif, y = proba)) + geom_boxplot()
classifFEM <- apply(resfunFEM$P, 1, which.max)
plotprofilsmoy(Velibdata, clustering = classifFEM, plotsem = FALSE)
qmplot(longitude, latitude, data = station, colour = as.factor(classifFEM))
stationcaract(station, classifFEM)
table(resICL$classification, classifFEM)
adjustedRandIndex(resICL$classification, classifFEM)
table(resICLbis$classification, classifFEM)
adjustedRandIndex(resICLbis$classification, classifFEM)